<html>
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
    <title>NDcost</title>
  </head>
  <body bgcolor="#FFFFFF">
    <center>Scilab Function</center>
    <div align="right">Last update : 14/03/2005</div>
    <p>
      <b>NDcost</b> -  generic external for optim
  computing gradient using finite differences</p>
    <h3>
      <font color="blue">Calling Sequence</font>
    </h3>
    <dl>
      <dd>
        <tt>[f,g,ind]=NDcost(x,ind,fun,varargin) </tt>
      </dd>
    </dl>
    <h3>
      <font color="blue">Parameters</font>
    </h3>
    <ul>
      <li>
        <tt>
          <b>x</b>
        </tt>: real  vector or matrix</li>
      <li>
        <tt>
          <b>ind</b>
        </tt>: integer parameter (see optim)</li>
      <li>
        <tt>
          <b>fun</b>
        </tt>: Scilab function with calling sequence <tt>
          <b>F=fun(x,varargin)</b>
        </tt>
           varargin may be use to pass  parameters <tt>
          <b>p1,...pn</b>
        </tt>
      </li>
      <li>
        <tt>
          <b>f</b>
        </tt>: criterion value  at point <tt>
          <b>x</b>
        </tt> (see optim)</li>
      <li>
        <tt>
          <b>g</b>
        </tt>: gradient value  at point <tt>
          <b>x</b>
        </tt> (see optim)</li>
    </ul>
    <h3>
      <font color="blue">Description</font>
    </h3>
    <p>
  This function can be used as an external for <tt>
        <b>optim</b>
      </tt> to minimize
   problem where gradient is too complicated to be programmed. only the
   function <tt>
        <b>fun</b>
      </tt> which computes the criterion is required.</p>
    <p>
   This function should be used as follow:
   <tt>
        <b>[f,xopt,gopt]=optim(list(NDcost,fun,p1,...pn),x0,...)</b>
      </tt>
    </p>
    <h3>
      <font color="blue">Examples</font>
    </h3>
    <pre>
// example #1 (a simple one)
//function to minimize
function f=rosenbrock(x,varargin)
   p=varargin(1)
   f=1+sum( p*(x(2:$)-x(1:$-1)^2)^2 + (1-x(2:$))^2)
endfunction

x0=[1;2;3;4];
[f,xopt,gopt]=optim(list(NDcost,rosenbrock,200),x0)

// example #2: This example (by Rainer von Seggern) shows a quick (*) way to
//             identify the parameters of a linear differential equation with 
//             the help of scilab.
//             The model is a simple damped (linear) oscillator:
//
//               x''(t) + c x'(t) + k x(t) = 0 ,
// 
// and we write it as a system of two differential equations of first
// order with y(1) = x, and y(2) = x':
//
//     dy1/dt = y(2)
//     dy2/dt = -c*y(2) -k*y(1).
//
// We suppose to have m measurements of x (that is y(1)) at different times 
// t_obs(1), ..., t_obs(m) called x_obs(1), ..., x_obs(m)  (in this example
// these measuresments will be simulated), and we want to identify the parameters
// c and k by minimizing the sum of squared errors between x_obs and y1(t_obs,p).
// 
// (*) This method is not the most efficient but it is easy to implement.
// 
function dy = DEQ(t,y,p)
   // The rhs of our first order differential equation system.
   c =p(1);k=p(2)
   dy=[y(2);-c*y(2)-k*y(1)]
endfunction

function y=uN(p, t, t0, y0)
   // Numerical solution obtained with ode. (In this linear case an exact analytic
   // solution can easily be found, but ode would also work for "any" system.)
   // Note: the ode output must be an approximation of the solution at
   //       times given in the vector t=[t(1),...,t($)]  
   y = ode(y0,t0,t,list(DEQ,p))
endfunction

function r = cost_func(p, t_obs, x_obs, t0, y0) 
   // This is the function to be minimized, that is the sum of the squared
   // errors between what gives the model and the measuments.
   sol = uN(p, t_obs, t0, y0)
   e = sol(1,:) - x_obs
   r = sum(e.*e) 
endfunction

// Data
y0 = [10;0]; t0 = 0; // Initial conditions y0 for initial time t0. 
T = 30;  // Final time for the measurements.

// Here we simulate experimental data, (from which the parameters
// should be identified).
pe = [0.2;3];  // Exact parameters
m = 80;  t_obs = linspace(t0+2,T,m); // Observation times
// Noise: each measurement is supposed to have a (gaussian) random error
// of mean 0 and std deviation proportional to the magnitude
// of the value (sigma*|x_exact(t_obs(i))|).
sigma = 0.1;  
y_exact = uN(pe, t_obs, t0, y0);
x_obs = y_exact(1,:) + grand(1,m,"nor",0, sigma).*abs(y_exact(1,:));

// Initial guess parameters
p0 = [0.5 ; 5];  

// The value of the cost function before optimization:
cost0 = cost_func(p0, t_obs, x_obs, t0, y0); 
mprintf("\n\r The value of the cost function before optimization = %g \n\r",...
                                                                       cost0)

// Solution with optim
[costopt,popt]=optim(list(NDcost,cost_func, t_obs, x_obs, t0, y0),p0,...
                                                       'ar',40,40,1e-3);

mprintf("\n\r The value of the cost function after optimization  = %g",costopt)
mprintf("\n\r The identified values of the parameters: c = %g, k = %g \n\r",...
                                                               popt(1),popt(2))

// A small plot:
t = linspace(0,T,400);
y = uN(popt, t, t0, y0);
clf();
plot2d(t',y(1,:)',style=5)
plot2d(t_obs',x_obs(1,:)',style=-5)
legend(["model","measurements"]);
xtitle("Least square fit to identify ode parameters")
 
  </pre>
    <h3>
      <font color="blue">See Also</font>
    </h3>
    <p>
      <a href="optim.htm">
        <tt>
          <b>optim</b>
        </tt>
      </a>,&nbsp;&nbsp;<a href="../programming/external.htm">
        <tt>
          <b>external</b>
        </tt>
      </a>,&nbsp;&nbsp;<a href="derivative.htm">
        <tt>
          <b>derivative</b>
        </tt>
      </a>,&nbsp;&nbsp;</p>
  </body>
</html>
